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We studied the role that electron gas has on the formation of nuclear structures at subsaturation densities and 
low temperatures (7' « 1 MeV). Using a classical molecular dynamics model we studied isospin symmetric and 
asymmetric matter at subsaturation densities and low temperatures varying the Coulomb interaction strength. 
The effect of such variation was quantified on the fragment size multiplicity, the inter-particle distance, the 
isospin content of the clusters, the nucleon mobility and cluster persistence, and on the nuclear structure shapes. 
We found that the presence of an electron gas distributes matter more evenly, disrupts the formation of larger 
objects, reduces the isospin content, and modifies the nucleon average displacement, but does not affect the 
inter-nucleon distance in clusters. The nuclear structures are also found to change shapes by different degrees 
depending on their isospin content, temperature and density. 
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I. INTRODUCTION 


Understanding nuclear matter at subcritical temperatures 
and subsaturation densities is of importance to understand 
phenomena ranging from heavy-ion collisions to the stability 
of stellar objects. This region has been investigated through 
nuclear reactions leading to the fragmentation of the nuclei, 
and it is now known that excited nuclear matter can exhibit a 
first-order phase transition [1], spinodal decomposition [2, 3], 
critical behavior [4] and other phenomenon associated with 
liquid-gas systems. More recently, such studies have been 
extended to neutron-rich nuclei [5-8]. 

At lower temperatures (T < 1MeV ) and subsatura- 
tion densities, theoretical studies predict that nuclear matter 
forms non-homogeneous structures known as nuclear pasta. 
Synoptically, static models [9—14] (e.g., liquid-drop model, 
Thomas-Fermi, and Hartree-Fock) find such structures from a 
balance between the different components of the nuclear and 
Coulomb energies, while dynamical models [15—24] (e.g., 
quantum and classical molecular dynamics) find them as e- 
quilibrium states of dynamical evolutions of such systems. 

The non-homogeneous phases are expected to exist on the 
outer crust of neutron stars and play a role in the emission 
of neutrinos [21], starquakes and pulsar glitches [25]. In 
such stellar conditions an unbalance between the weak de- 
cay n — p+e c v, and its opposite process, p +e —^ n +v 
is expected to yield an accumulation of neutrons and pro- 
duce a neutralizing electron gas. Thus researchers have s- 
tudied such nuclear structures while embedded in an electron 
gas. Although most calculations have shown that the struc- 
tures are formed thanks to the interplay of nuclear forces with 
the Coulomb forces of the electron gas, none have dissected 
this effect quantitatively. The question we aim to answer is if 
the nuclear pasta can structures form without the embedding 
electron gas? 
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This investigation will focus on the role of the Coulomb 
forces on the creation of the pasta using a classical molecular 
dynamics which was shown to be effective in these type of 
studies [22]. By using a microscopic view of the nucleons 
and a combination of topological tools the different struc- 
tures formed can be characterized into recognizable patterns, 
opening the door to a quantifiable analysis of the effect the 
Coulomb forces of the electron gas have on the pasta forma- 
tion. 

In 2003 [24] a static liquid-drop model was used to study 
the effect of an electron gas on cold nuclear structures, and 
it was found that its main effect was to extend the range of 
densities where bubbles and clusters appear. Later, in a 2005 
study [26] a density functional method was used (at zero tem- 
perature) to find that the density region in which the pasta 
exists becomes broader when the electron gas is taken into 
account. The authors remarked the importance of extending 
such study to finite temperatures and with dynamical models. 

Recently, Monte Carlo simulations, which approximated 
the Coulomb interaction via an Ewald summation [21, 27— 
29], determined the lowest proton fractions that are compati- 
ble with 8 equilibrium in neutron star crust environments (at 
T = 1.0 MeV). The authors recommended that further inves- 
tigation be performed with molecular dynamics simulations 
to study dynamical observables and to examine the role of 
the temperature in the study of nuclear matter under subsatu- 
ration density and low temperature conditions. 

The goal of this work is to study the effect that electron gas 
has on the formation of the nuclear pasta using a dynamical 
model, at non-zero temperatures, and in a systematic way 
that takes advantage of the tools developed in a previous 
work [22]. It must be remarked that it is not the intention 
of this work to study either the structure or other properties of 
neutron star crusts nor to propose the present molecular dy- 
namics method for such purpose; the only goal of the present 
exercise is —again- to understand the role the electron gas 
plays in determining the pasta shapes. This will be achieved 
through the calculation of the pasta structures with and with- 
out the effect of the electron gas as well as with softened 
interactions which, of course, do not exist in nature. Reit- 
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Fig. 1. (Color online) Characteristic structures obtained with (a) and without (b) the proton-electron gas interaction for symmetric matter 


(x = 0.5) at density p = 0.015 fm? and temperature T = 0.1MeV. 


erating, we are not hoping to mimic neutron star crusts nor 
proposing the existence of nuclear structures with weakened 
electron gas interactions. 

In the next section, the model used (classical molecular 
dynamics) is briefly introduced along with the analysis tools 
employed. Section III presents a series of results of the sim- 
ulations and analyses performed. We summarize some of our 
findings in Section IV and present a quick look at work in 
progress. 


Il. MODEL USED 


The model used here was developed to study nuclear reac- 
tions from a microscopic point of view, and it is composed 
of a molecular dynamics code retrofitted with cluster recog- 
nition algorithms and a plethora of analysis tools. The model 
was first used for infinite nuclear matter elsewhere [8, 22, 23], 
here we simply mention some basic ingredients of the model. 

The classical molecular dynamics model (CMD), first in- 
troduced in [30], has been used in heavy-ion reaction studies 
to explain experimental data [31], identify phase transitions 
and other critical phenomena [32-35], explore the caloric 
curve of nuclear matter [36, 37] and isoscaling [38, 39]. CMD 
uses two two-body potentials to describe the motion of “nu- 
cleons" by solving their classical equations of motion. The 
potentials were developed phenomenologically by Pandhari- 
pande [40-42], we direct interested readers to this original 
work for a review of the basic properties of nuclear matter 
and finite nuclei that can be described by this model. The 
potential used is 


Vap(r) =V; [exp(—prr)/r — exp(—Hrre)/re] 
— Va [exp(—Har)/r — exp(—Harc)/re], (1) 
Vyn (r) = Vo lexp(—uor)/r — exp(—Hore)/re] , 


where Vj, is the potential between a neutron and a pro- 
ton, and Vy vy is the repulsive interaction between either nn 
or pp. The cutoff radius is re = 5.4fm and for r > re 
both potentials are set to zero. The Yukawa parameters pr, 
Ha and uo were determined by Pandharipande to yield an 
equilibrium density of pọ = 0.16fm ?, a binding energy 
E(po) = 16 MeV/nucleon and a selectable compressibility of 
250 MeV for the “Medium” model, and 535 MeV for the “S- 
tiff” [40-42]. We remark that, since the potential Vy y does 
not permit bound states between identical nucleons, pure neu- 
tron matter is unbound and, differing from potentials used by 
other models [21], the Pandharipande potentials have a hard 
core. 

The main advantage of the CMD model is the possibility 
of knowing the position and momentum of all particles at all 
times, as this allows the study of the structure of the nuclear 
medium from a microscopic point of view. The outcome of 
CMD, namely, the time evolution of the particles in (r.p) can 
be used to study the resulting pasta structures. 

On the other hand, the CMD model fails to comply with 
several quantum effects. Luckily, the model can still be 
applied for the nucleon dynamics since, as explained else- 
where [22, 43, 44], the mean inter-particle distance (V/N)!/? 
in nuclear matter in the region of interest (p < 0.1fm~°, 
T < 1 MeV) is larger than the nucleon thermal wave- 
length [8]. Other flaws of the CMD are to be discussed after 
the results are presented. 

To simulate an infinite medium, systems with thousands 
of nucleons were constructed using CMD under periodic 
boundary conditions. Cases symmetric in isospin (i.e. with 
x = Z/A = 0.5, 1000 protons and 1000 neutrons), and 
cases with asymmetric isospin (with x = 0.3 comprised of 
1000 protons and 2000 neutrons) were constructed in cu- 
bical boxes with sizes adjusted to have densities between 
0.01fm^? < p < po; it must be remarked that by fixing 
the mass and size of each cell, the scale of the density fluctu- 
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ations is also determined. Our previous study [22], presented 
a smórgásboard of structures obtained through this method. 


A. Coulomb interaction 


To study nuclear matter embedded in a degenerate elec- 
tron gas, the system is first taken as overall neutral (8- 
equilibrated), and then a proper approximation is used to take 
into account the Coulomb interaction (which is of infinite 
range). The two most common approaches are the Thomas- 
Fermi screened Coulomb potential (used, e.g., in QMD [16]) 
and the Ewald summation procedure [45]; CMD has been 
used with both approximations (see [22] for a comparison of 
methods under CMD) but for computational reasons we opt 
for the use of a screened Coulomb potential. 

We approximate the electron gas as a uniform ideal Fermi 
gas with the same number density as the protons. Its inclu- 
sion in the CMD is through the use of a screened Coulomb 


potential obtained from the Poisson equation, VE) (r) = 
(e?/r) exp(—r/A), where the screening length, A, is related 
to the electron mass, Fermi momentum, and number densi- 
ty. Here, pragmatically, and after testing various lengths, we 
used the prescription used before [16, 21] and set A=10 fm 
which satisfies the requirement of being sufficiently smaller 


than the size of the simulation cell, L — (A/ p) s 


B. Simulation procedure 


The trajectories of the nucleons are then governed by the 
Pandharipande and the screened Coulomb potentials. The 
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Fig. 2. Effect of the screened potential on the fragment size multi- 
plicity. Plots show the distribution of cluster sizes observed in 200 
configurations of x = 0.3 nuclear matter at density p = 0.015 fm^? 
and temperatures T = 1.0 MeV (top) and T = 0.1 MeV (bot- 
tom). The figure on the left panel correspond to configurations with 
Coulomb interaction, and those on the right to the case without such 
potential. 
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X-0.5, p = 0.072 fm, T=0.1 MeV 


Fig. 3. Examples of the radial correlation function for varying 
strengths of the Coulomb potential: œ = 1 (full Coulomb), 0.8, 
0.2, and 0 (without Coulomb). 


equations of motion are solved using a symplectic Verlet inte- 
gration with energy conservation of O(0.01%). The nuclear 
system is force-heated or cooled using isothermal molecular 
dynamics with the Andersen thermostat procedure [46] with 
small temperature steps. We are interested in the range of 
0.1 MeV < T < 1.0MeV and densities 0.01 fm^? < ps 
po. In comparison to heavy-ion reactions, these conditions 
correspond to practically semi-frozen states. 

To obtain reliable statistics we sample each T', p and x con- 
figuration 200 times in a stationary ergodic process. That is, 
after reaching equilibrium at a given set of values of T, p and 
x, the (r, p) values of all nucleons are stored for further pro- 
cessing while the system continues its evolution for a “long” 
time (4000 time steps) until it reaches a state independent of 
the previous one, point at which the procedure is repeated 
for a total of 200 times. The recorded nucleon positions and 
momenta are then used later to identify clusters and to char- 
acterize the structure as explained next. 

Given that, for the densities and temperatures of interest, 
the nucleons do not have large mobilities, the task of identify- 
ing self-bound cumuli is sufficiently simple as to be achieved 
with a basic "Minimum Spanning Tree" (MST) algorithm. 
The MST finds cluster membership simply by identifying 
those particles that are closer to each other than some clusteri- 
zation radius, which here is set to 3.0fm based on the range of 
the nuclear potential and densities investigated. The generic 
method was modified to recognize fragments that extend into 
any of the adjacent periodic cells. 


C. Analysis tools 


Once the fragments have been identified, their global prop- 
erties are studied as follows. MST yields information about 
the fragment multiplicity. In Ref. [22] systems at x — 0.3 and 
0.5 and at various densities (0.01 fm^? < p < 0.015 fm-?) 
and temperatures (0.1 MeV < T < 1.0MeV) exhibited 
increases in multiplicity with temperature (in opposition to 
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Fig. 4. x content of the clusters formed in 200 configurations of 
asymmetric matter (x = 0.3) at density p = 0.015 fm ? and tem- 
perature 7' = 1.0 MeV. The top panel shows the case with Coulomb 
and the bottom one the case without such interaction. Notice that the 
abscissas have different scales. 


the rise-fall that takes place in intermediate-energy heavy-ion 
reactions), as well as the formation of percolating structures 
that encompass most of the nucleons in the simulation. 


At a microscopic level, the dynamics of the nucleons in 
such systems can be quantified through their average dis- 
placement as a function of “time”, i.e., through the time steps 
of the simulation; for the cases mentioned before [22], the 
mobility between configurations was found to be between 2 
fm and 20 fm. Likewise, the microscopic stability of the clus- 
ters can be gauged through the “persistency” [43] which mea- 
sures the tendency of members of a given cluster to remain 
in the same cluster. For the cases previously studied, cold- 
er clusters (T' < 0.5 MeV) exchange less than half of their 
nucleons during 200 ergodic configurations, while those at 
higher temperatures (T > 0.8 MeV) replace up to 90% of 
their constituents. 


Another interesting descriptor is the isospin content of each 
cluster produced, i.e., their z values. In the previous study it 
was found that in the case of x = 0.3, small clusters (A < 
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Fig. 5. RM S displacement of nucleons during the evolution of the 
stationary ergodic process with Coulomb (continuous line) and with- 
out Coulomb (dashed line) as labeled; see text for details. Notice that 
at low temperatures the mobility with Coulomb is smaller than with- 
out Coulomb, while at higher temperatures the effect is the opposite. 


10) tend to have less protons than the global average, while 
such effect was not present in the case of x = 0.5. Another 
global characterization tool is the pair correlation function, 
g(r), which gives information about the spatial ordering of 
the nuclear medium by comparing the average local density 
to the global density. In the previous study g(r) showed that 
nucleons in clusters have an interparticle distance of about 
1.7fm at all studied densities, and was useful in determining 
the onset of lasagna-like structures. 


Proceeding beyond global measures, the shapes of nuclear 
structures can be characterized by their volume, surface area, 
mean curvature, plus an interesting construct known as the 
Euler characteristic. These four objects comprise what is 
known as the “Minkowsky functionals” which completely 
describe all topological properties of any three-dimensional 
object. The evaluation of the mean curvature (averaged over 
the whole surface of the object) and Euler number (a topo- 
logical invariant connected to the vertices, edges and faces 
of polyhedra), however, requires the mapping of the nucle- 
ar clusters into a polyhedra. This procedure can be accom- 
plished through the Michielsen—De Raed algorithm [19, 47]. 


In [22], it was shown that generic structures, such as “gnoc- 
chi” » “spaghetti”, “lasagna” and “crossed-lasagnas” or “jun- 
gle gym” and their inverse structures (with voids replacing 
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Fig. 6. Persistence with and without Coulomb as it evolves through 
200 configurations of the stationary ergodic process. 


particles and vice versa), all have well defined and distinct 
values of the mean curvature and Euler numbers. It should be 
remarked that the assignment of curvature and Euler values to 
a given structure becomes relative to the overall scale. This 
becomes important for structures with near-zero Euler num- 
ber which signal spaghetti, lasagna and their anti-structures. 
Curled spaghetti rods of a given length, for instance, can eas- 
ily oscillate between positive and negative near zero Euler 
numbers by curling or uncurling a segment. This classifica- 
tion, which was respected at all x values, densities and tem- 
peratures studied, is shown in Table 1. 

It is with this model and arsenal of tools that we now 
embark on a study of the effect that electron gas has on the 
formation of the nuclear pasta. 


TABLE 1. Classification Curvature - Euler 
Curvature < 0 Curvature ~ 0 Curvature > 0 


Euler > 0 Anti-Gnocchi Gnocchi 
Euler ~ 0 Anti-Spaghetti Lasagna Spaghetti 
Euler < 0 Anti-Jungle Gym Jungle Gym 


IIl. THE EFFECT OF THE ELECTRON GAS 


Along the lines of previous studies [24, 26, 29] and adher- 
ing to their advise, in this study we used a molecular dynam- 
ics model to study spatial and dynamical effects the electron 
gas may produce on the nuclear structures at non-zero temper- 
atures. Specifically we extended the approach of one of the 
previous studies and compared pasta structures obtained with 
and without the electron gas as well as with varying strengths. 

In summary, pasta structures were obtained at subsatura- 
tion densities (0.015fm^? < p < 0.072fm-?) and low 
temperatures (0.1 MeV < T' € 1.0MeV) for symmetric 
(x = 0.5) and asymmetric (x = 0.3) cases. Altogether 
200 simulations were carried out per each combination of 
(o, T, x). To see the effect of the electron gas, each case was 
"cooked" repeatedly with a screened Coulomb potential with 
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a varying amplitude, i.e. with aV em. where 0 € o < 1. 


For a fair comparison, all corresponding cases with different 
values of a were produced with identical initial conditions of 
{x, p, T) and started off from the same initial random config- 
uration. As an illustration, Fig. 1 shows two corresponding 
structures obtained with and without the electron gas. 

The effect of the screened potential can also be observed on 
the fragment size multiplicity. Fig. 2 shows the distribution 
of cluster sizes observed in 200 configurations of asymmetric 
matter (x = 0.3) at density p = 0.015 fm? and temperatures 
T = 0.1MeV and T' = 1.0MeV. The figures correspond 
to configurations with and without the Coulomb interaction. 
The increase in the number of fragments of sizes A > 500 as 
a goes from a = 1 to a = 0 underline the role the electron 
gas has on determining the mass distribution. 

The change of the inner structure can be quantified through 
the use of the radial distribution function, g(r). An example 
of this effect is shown in Fig. 3, which shows this function cal- 
culated with the structures obtained with symmetric nuclear 
matter (x = 0.5) at density p = 0.072 fm-? and tempera- 
ture T' — 0.1 MeV and for varying strengths of the Coulomb 
potential, namely o — 1 (full Coulomb), 0.8, 0.2, and 0 (with- 
out Coulomb). As it can be easily observed, the electron gas 
does not affect the nearest neighbor distances, although the 
magnitude of the function for the second nearest neighbors 
increases as Coulomb weakens. 

The isotopic content on the cluster is also affected by the 
presence of electron gas. Fig. 4 shows the x content of the 
clusters formed in 200 configurations of asymmetric matter 
(x = 0.3) at density p = 0.015 fm? and temperature T = 
1.0MeV. Clearly visible is the enhancement up to x ~% 0.4 
that would occur without Coulomb. Also noticeable is the 
increase of the maximum fragment sizes which grow up to 
A ~ 800 without the electric interaction. 

As suspected in [26], in this study we found that the 
interactions between protons and electrons are responsible 
for the rearrangement of protons, although we found this to 
be of varying degree. Fig. 5 shows the RMS displacement 
of nucleons as a function of the evolution of the stationary 
ergodic process. Each curve is made of 200 points, each of 
which represents the average displacement (with respect to 
the original configuration) of the 2000 nucleons used in the 
x = 0.5 cases shown. For each case tracked (p = 0.015 
and 0.72 fm-? at T = 0.1 and 1.0 MeV), the displacements 
obtained with and without Coulomb are presented. Curiously, 
at low temperatures Coulomb seems to enhance the nucleon 
mobility much more than at high temperatures. 

A similar effect is observed in the persistence, which is the 
tendency of nucleons in a given fragment to remain in the 
same cluster. Fig. 6 shows the persistence with and without 
Coulomb as it evolves through 200 configurations of the sta- 
tionary ergodic process, measured with respect to the original 
configuration of 3000 nucleons with z — 0.3, p — 0.015 and 
T = 1.0MeV. As seen in the mobility, Coulomb enhances 
the transfer of nucleons thus decreasing the persistence as 
shown by the two curves. 

The effect of electron gas on the shape of the structures 
formed can be quantified using the Euler Characteristic and 
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Mean Curvature introduced before. Following the procedure 
explained in more detail elsewhere [22], digitized polyhe- 
dra were constructed for each of the 200 nuclear structures 
obtained with and without Coulomb at each value of (x, p, T') 
considered. The corresponding values of the mean curvature 
and Euler numbers were calculated and, although the config- 
urations corresponding to the same set of conditions were not 
identical, their curvature and Euler numbers showed robust 
average values. 

The actual values of the curvature in all the studied configu- 
rations varied approximately from —6,000 to 12,000 and those 
of the Euler numbers from —700 to 1,000. As these values 
depend on the characteristics of the polyhedra (i.e., number 
of sides, corners, etc.) they yield only relative information. 
Table 2 shows the curvature-Euler numbers as percentages 
of the maximum values obtained, and Fig. 7 shows their loca- 
tion on the curvature-Euler plane as well as their changes in 
position as the Coulomb strength is diminished. The standard 
deviations of the points are, roughly speaking, of the size of 
the points used in the plots. 


TABLE 2. Classification Curvature - Euler 


p T a=1 a=0 

(fm^?) (MeV) Curvature Euler Curvature Euler 
A 0.072 1.0 -14.5 8.9 -0.2 0.47 
B 0.072 0.1 -3.18 -1.58 0.245 -0.19 
C 0.015 0.1 14.76 0.89 4.28 0.35 
D 0.015 1.0 18.43 -0.17 2.28 0.86 
E 0.015 1.0 75.75 -61.2 91.76 -24.3 
F 0.015 0.1 90.03 45.4 95.3 15.9 
G 0.072 1.0 -24.74 44.7 -37.6 60.2 
H 0.072 0.1 -37.38 79.8 -49.28 93.3 

A. Analysis 


The results obtained in this study allow us to understand 
the effect that electron gas has on the formation of the nuclear 
pasta. Perhaps the most interesting result of this investigation 
is the fact that the characteristic topological structures of the 
nuclear pasta, which are usually associated with the presence 
of competing interactions of the long range Coulomb poten- 
tial and short range nuclear force- exist even when there is 
no Coulomb (i.e., when a — 0). We believe this to be the 
result of the competition between attractive and repulsive nu- 
clear forces at very low temperatures. In our case this is due 
to the attraction of the V,,, and the repulsion of the Vpn and 
Vpp which operates at medium ranges. Incidentally, we find 
the same effect also in other potentials [29] as was reported 
elsewhere [23]. 

The role of the Coulomb forces introduced by the elec- 
tron gas is to emphasize such competition between attractive- 
repulsive forces. The examples of structures formed with and 
without Coulomb (Fig. 1) suggest that the electric forces tend 
to distribute matter more and form less-compact objects. This 
behavior is in agreement with the reduction of the fragment 
size multiplicity as illustrated in Figs. 2 and 4. 
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Fig. 7. Average values of the Curvature and Euler numbers of 
the structures listed in Table 2; circles correspond to structures 
with Coulomb and squares to structures without Coulomb, arrows 
indicate the average displacement of the structures as o goes from 1 
to 0. 


At a microscopic scale, Fig. 3 shows that varying the 
strength of the Coulomb interaction does not change the inter- 
particle distance, but certainly decreases a bit the x content 
of the fragments as can be seen in Fig. 4. In agreement with 
this, the increment of nucleon mobility produced by Coulomb 
(Fig. 5) gets reflected in a reduction of the persistence (Fig. 6). 


With respect to the modification of the shapes due to the 
electric interaction, we find that Coulomb tends to reduce 
both the curvature and Euler number of the structures. Indeed, 
Coulomb appears to “compactify” symmetric matter into 
spaghetti and lasagna type structures as shown clearly in 
Fig. 7 for symmetric matter (top panel). Of course, the effect 
becomes blurrier in systems with smaller charge densities, 
such as the asymmetric matter case of x = 0.3 (bottom pan- 
el). 

After having seen the results of the present study, it is now 
easier to gauge the effect of those quantum features not in- 
cluded in our classical model. CMD is at best an independent- 
particle approximation to a quantum system in which all high- 
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er order correlations are ignored. This excludes all effects 
ascribed to the closing of shells (magic numbers, shell po- 
larization and deformations) and to any collective phenom- 
ena (pairing rotations and vibrations and their effect on fis- 
sion decays, formation of Cooper pairs and their superfluidity 
effects). Although some of these effects are expected to be 
noticeable in the low density low temperature region of our s- 
tudy, it is easy to see that few of them —if any- have a bearing 
on the macroscopic structure of the clusters formed. Like- 
wise, as all of the missing quantum effects are independent of 
the existence (or lack) of an embedding electron gas, it stands 
to reason that they would not have any effect on the result- 
s. In other words, had the comparison between pasta shapes 
with and without Coulomb interactions been performed with 
a model that included all quantum effects, the results would 
have been the same. 


IV. CONCLUDING REMARKS 


The effect of electron gas on symmetric and neutron-rich 
matter was studied at low densities and temperatures by vary- 
ing the Coulomb interaction strength. Its effects on the frag- 
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